module module_slab
use polymorphic_complextaylor
implicit none
integer nx,i,no,nv
type(real_8) y,f,n,x,y0,fac
real(dp) xr,yr,dyr,dtheta,fc,radius_of_curvature,fr
integer,  parameter  :: npo=200
integer :: nob=2
type dpoint
 real(sp) x,y
end type dpoint
type(dpoint) point_surface(-npo:npo)
type(dpoint) intercept_axis(-npo:npo)

contains

subroutine compute_surface
implicit none

open(unit=10,file='lens_shape.txt')
open(unit=11,file='intercept.txt')

fac=0.85d0

no=nob

nv=1;
call init(no,nv)  
call alloc(y,f,n,x,y0)
fr=0.9d0
f=fr
n=1.5d0

call integrate_y(y,f,n)
radius_of_curvature=1.d0/(-2.d0*(y.sub.'2'))
!write(6,*) " radius of curvature ", radius_of_curvature

xr=fac*radius_of_curvature/npo
do i=-npo,npo
x=i*xr
call compute_position_y(y,x,y0)
 yr=y0
 point_surface(i)%y=i*xr
 point_surface(i)%x=yr

!write(6,*) " yr ", yr
enddo



xr=fac*radius_of_curvature/npo
do i=-npo,npo
if(i/=0) then
x=i*xr
call compute_position_yp(y,x,y0)
 dyr=y0
call compute_position_y(y,x,y0)
 yr=y0
 dtheta= asin(dyr/sqrt(1.d0+dyr**2))-asin(n*dyr/sqrt(1.d0+dyr**2))
 fc=(i*xr+yr*tan(dtheta))/tan(dtheta)
 intercept_axis(i)%y=0.d0
 intercept_axis(i)%x=fc
 else
 intercept_axis(i)%y=0.d0
 intercept_axis(i)%x=f
 endif
 !write(6,*) fc
enddo


call kill(y,f,n,x,y0)

!goto 111

do i=-npo,npo

write(10,*) point_surface(i)%x,point_surface(i)%y
write(11,*) intercept_axis(i)%x,intercept_axis(i)%y

enddo

close(10);close(11);
end subroutine compute_surface

subroutine integrate_y(y,f,n)
type(real_8) y,yb,f,x,n,dydx
integer order,j
integer, allocatable :: ju(:)

allocate(ju(c_%nv))
ju=0

call alloc(x,dydx,yb)

yb=0.d0

x=morph(1.d0.mono.'1')
! first compute derivative

do i=1,c_%no
dydx=x/(  (f-yb)-n*sqrt((f-yb)**2+x**2) )

y=0.d0
do j=0,c_%no
  ju(1)=j
  y=y+ (dydx.sub.ju)*x**(ju(1)+1)/(ju(1)+1)
enddo
yb=y
!call print(yb,6)
enddo


deallocate(ju)
call kill(x,dydx)
end subroutine integrate_y

subroutine compute_position_y(y,x,y0)
type(real_8) y,x,y0
integer i
integer, allocatable :: ju(:)

allocate(ju(c_%nv))
ju=0

y0=0.d0

do i=0,c_%no
ju(1)=i
y0=y0+ (y.sub.ju)*x**i
enddo


deallocate(ju)

end subroutine compute_position_y

subroutine compute_position_yp(y,x,y0)
type(real_8) y,x,y0,dydx
integer i
integer, allocatable :: ju(:)

call alloc(dydx)
allocate(ju(c_%nv))

dydx=y%t.d.1

ju=0

y0=0.d0

do i=0,c_%no
ju(1)=i
y0=y0+ (dydx.sub.ju)*x**i
enddo


deallocate(ju)
call kill(dydx)

end subroutine compute_position_yp



end module module_slab
